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PROJECT SUMMARY 

Radiation heat transfer in rocket motor combustion chambers and plumes has 
historically been analyzed by simulating the flowfield without considering radiation, 
and as a second step calculating the radiation without modifying the flowfield 
solution. New and advanced propulsion system concepts have made these thermal 
fields more severe, so that the uncoupled solutions are no longer acceptable. This 
study coupled the radiative transport directly to the conservation equations and, 
thereby, yielded a simultaneous solution to all of the governing equations. 

The three commonly used propellant systems: H 2 /0 2 , RP-1/0 2 , and solid 
propellants radiate as: molecular emitters, non-scattering small particles, and 
scattering larger particles, respectively. This study addressed all of these types of 
radiation. Methods of representing the spectral character of these radiators were 
reviewed to select the most efficient radiation models for use in a coupled solution. 
Line-of-sight (LOS) calculations were made to verify that the radiation models selected 
were sufficiently accurate for radiation heat transfer predictions. 

The nature of the radiation/convection coupled problem is that the divergence 
of the radiative heat flux must be included in the energy equation and that the local, 
volume averaged intensity must be determined by a solution to the radiative transfer 
equation (RTE). For a rigorous solution, the RTE must be solved along several LOS 
at each point in the flowfield at each computation step. Such a procedure is 
extremely costly; therefore, approximations are needed. The ordinary and improved 
differential approximations were developed to provide the required solution. Both of 
these solutions were provided as additions to the FDNS code (a Navier-Stokes solver 
for turbulent, two-phase combusting flow) to form the radiation/convection coupled 
solution. The ordinary differential approximation is applicable to combustion chambers 
where the flowfield is optically dense, and the improved differential approximation is 
applicable to nozzle flow and plumes which are of intermediate optical depth. 

The thermal analysis provided by this radiation/convection coupled methodology 
is applicable to rocket motor and vehicle base heating design. It is of particular value 
in analyzing ablator material response in solid rocket motors and in gas/solid 
combustion phenomena in hybrid rocket motors, where the radiation heat load is 
critical to the system performance. Industrial applications to coal-fired furnace and 
glass melting-pot analyses are also feasible. 
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RADIATION/CONVECTION COUPLING 
IN ROCKET MOTOR & PLUME ANALYSIS 

1 .0 Introduction 

A computational fluid dynamics (CFD) code is presented which fully couples 
the radiation/convection interaction. Radiation heat flux is included in the energy 
equation, and the radiation transport equation is solved. These radiative terms and 
equations were added to the FDNS CFD code. Justification for choosing the 
radiation methods presented herein are given, along with the limitations of these 
methods, for various motor types. The code is flexible in that higher order 
radiative approximations can be added to the current order, as needed. 

1.1 Methodology 

Radiation within a flowfield depends upon the absorption and emission of 
electromagnetic energy to or from gaseous and particulate species. In the case of 
particles, there is an additional effect caused by the scattering of photons toward 
or away from a point, from all locations and directions within the flowfield. To 
account for these effects, the first step was to determine the radiative properties 
of the constituent species, namely the absorption and scattering cross-sections. 
The next step required the solution to an approximation to the integro-differential 
radiative transport equation (RTE). This resulted in two approximate methods, one 
valid for optically thick regions and the other valid for arbitrary optical thicknesses. 
Finally, the radiative source term is determined and added to the right-hand-side of 
the energy equation, providing the coupling between the radiation and CFD solvers. 
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2.0 Absorption and Scattering Coefficients 

The three types of rocket motors currently being used - H 2 /0 2 , RP-1/0 2 , and 
solid motors - radiate in the following manner: 

1. Molecular radiation - infrared, non-luminous gaseous radiation associated 
mainly with rotational and vibrational energy modes; important for SSME 
nozzles with H 2 0 products. 

2. Non-scattering small particles, comprised mostly of soot, and applicable to 
RP-1/0 2 engines. The molecular radiators: C0 2 and H 2 0 are also important 
in these engines. 

3. Scattering large particles; specifically, those associated with solid 
propellants with Al 2 0 3 products. Again, molecular radiators are also 
important. 

This project addressed the radiation associated with H 2 0, C0 2 and Al 2 0 3 . 
Absorption and emission by soot will be added at a later date. 

2.1 Molecular Radiation 

Two methods were coded for the gaseous radiators - the total emissivity or 
mean band model and the wide band model. The former method assumes that the 
spectral width of each band is narrow, allowing an essentially constant Planck 
function over the range of each band so that all the bands are integrated to 
evaluate a total emissivity. Although band overlap is considered negligible, 
provisions are made to correct for this effect. However, since the emissivity is 
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obtained as an average over all bands, this method was not used for the gas 
radiation properties. The wide band model is a more rigorous procedure for finding 
absorption coefficients, and it also provides these absorption coefficients as a 
function of wavelength. This is an important requirement since the partial 
differential equation approximation obtained from the RTE is solved for a single 
wavelength. In the wide band model the entire spectral range of each band is 
modeled. In particular the exponential wide band model was coded. Band 
strengths and widths are related to parameters obtained from statistical mechanics. 
The correlation parameters used for H 2 0 and C0 2 (Refs. 1 and 2) are shown in 
Table 1. Five wavelengths for H 2 0 and six for C0 2 were deemed important for the 
analysis: 


^H20 — 71 , 6.3, 2.7, 1 .87, and 1 .38 jjm 
^co 2 =15, 10.4, 9.4, 4.3, 2.7 and 2.0 //m 

Narrow band models, such as those used in the JANNAF standard plume 
radiation code, SIRRM (Ref. 3) are the most accurate molecular radiation models 
currently available. However, the spectral integration required to use such models 
is too computationally intensive for use in a coupled solution. A comparison of 
total radiation for the SSME was made for the mean, wide, and narrow band (from 
SIRRM computations) models, and is shown in Table 2. A line-of-sight calculation 
was used to determine the SIRRM results. The emissivity computed from the 
mean band model is clearly inferior to the other two spectral methods. Differences 
between the wide and narrow band models are evident at the higher pressures, but 
the experimental molecular band model data in SIRRM were not collected at these 
high pressures. Since some of the wide band model data were taken at higher 
pressures, this method of correlation was chosen for further study. Comparison of 
these methods to a blackbody confirms the spectral nature of radiation throughout 
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Table 1 


Wide band model correlation parameters for various gases 
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the rocket motor. A similar comparison is made for gaseous products for the 
ASRM 48-5, in Table 3. 

2.2 Scattering (large) Particles 

Optical properties for Al 2 0 3 were obtained from Grumman's OPTROCK data 
(Ref. 4). Plots of the refractive index N, and the absorptive index N 2 are given in 
Figs. 1 and 2. These data were transformed to radiation properties by a Mie 
scattering code (Ref. 5) which determines absorption and scattering cross-sections 
as well as phase functions. The OPTROCK data includes wavelengths ranging 
from 0.2 to 25 /vm. Additional work has been reported (Ref. 6) to provide curve- 
fits of the absorptive index for the Al 2 0 3 optical property data. 

Scattering angle dependence on the diffractive, reflective and absorptive 
processes are determined through the phase function, P, as illustrated in Fig. 3. 
Note that a specific value of log 10 P is computed for each scattering angle 
COS (0). These figures indicate that as the size parameter (X = nD/A) increases, 
forward scattering increases. It is known that Mie scattering is a good 
approximation for forward scattering for large non-spherical particles. Since 
extinction is dominated by scattering in the forward direction (as seen in the 
figures), then scattering is not very sensitive to particle shape. This justifies using 
the Mie theory. Comparison of radiation properties between SIRRM and the Mie 
code is shown in Table 4. 

A line-of-sight (LOS) code was also used to determine radiance from 
particles. Comparison between SECA's two-flux model (Ref. 7) and the two-flux 
method in SIRRM appear in Fig. 4 for SSME and ASRM nozzles. A wavelength 
range from 0.6 to 8.0 //m was used in the two-flux method. The results indicate 
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that SECA's two-flux method is essentially identical to that in SIRRM. Although 
the two-flux method could be used for the present coupling analysis, many LOS 
calculations would be required to evaluate the volume-averaged intensity at each 
grid point. To reduce the computations required, the differential approximation to 
the RTE was chosen. 


3.0 Radiative Transport Equation 


The nature of the radiation/convection coupled problem is that the 
divergence of the radiative heat flux, "q r , must be included in the energy equation 
and that the local spectral intensity, I,,, of the radiation must be determined by a 
solution to the RTE (Ref. 8). The divergence of the radiative flux is given by: 


vg r = [4 7i I Xb ( T ) -j I x (r, fi) dCl'} dk 

O 

The RTE is given by: 

(f, G) = K x I Xb (T) - (k a + a x )I x (f,d) 


+ 


471 


[ P {i , Cl' -Cl) I x {r , Cl') dQ' 

J 4 it 


(l) 


( 2 ) 


where k and o are the absorption and scattering coefficients, respectively, l Ab is the 
Planck blackbody function, A is the wavelength, and r and S indicate position 

and direction, respectively, integrated over 4 n solid angles dCl ' . The l A is 
determined by solving Eq. 2 along several lines-of-sight for each point in the 
flowfield for each wavelength, or by using a further approximation. 
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4.0 Approximations to the RTE 
4.1 Ordinary Differential Approximation (ODA) 

The radiative intensity appears within the scattering integral of Eq. 2, 
resulting in an integro-differential equation. This results in a very costly system to 
solve. Consequently, the RTE was transformed to a simpler equation via a 
differential approximation. In this method the integral terms in Eq. 1 and Eq. 2 are 
evaluated by using the method of moments or spherical harmonics. Both the 
phase function and radiative intensity are rewritten as a generalized Fourier Series 
using Legendre Polynomials, with N terms in the spherical-harmonics expansion 
(the P N - approximation). The ordinary differential approximation (ODA) (Ref. 2, 9 
and 10), also known as the Pt - approximation, involves the solution to a partial 
differential equation to evaluate the integrated form of the RTE. Only the terms up 
to first order are retained in the series, resulting in: 

P{r, S' - fi) = 1 + Pl r-ca'-fi) <31 


V 


[1 -pju/3 


\ 

v G 

/ 


A 2 (4x/,-(3), 


(4) 


q. = — - — vG 
' 3 Yl <V, 


(5) 


where 

A 2 = 3 y 0 Ki 
Y 0 = 1 -Wj5 0 
Y\ = 1-wp,/3 
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( 7 ext is the extinction coefficient (absorption + scattering), ut is the albedo defined 
as the scattering coefficient divided by the extinction coefficient, p 0 and p, are the 
zero and first order scattering phase function parameters used in the expansion, 
respectively, and G is the incident radiation at a point 

G(r) = f l{TA')dCl l ( 6 ) 

For p 0 = 1 and p, = 0 the phase function defines linear isotropic scattering, while 
for p 0 and p, different from zero, linear anisotropic scattering is represented. Only 
linear isotropic scattering is coded at this time. 


i > 


A loss of generality results in that only optically thick media can be 
accurately described by the Pt - approximation. Since the combustor region of a 
rocket motor is optically thick, the Pt approximation was applied to that region. 
Writing Eq. 4 in cylindrical coordinates results in 


1 a ( r dG 


r dr 


dr 


1 + <?G 

r 2 5<t> 2 dr 2 


(7) 


-A 2 G = -4rc A 2 I b (T) 

For the axisymmetric nozzle, the azimuthal term (0) drops out. Although Eqs. 6 
and 7 are a function of wavelength, the A was removed for brevity. 


The P, - approximation eliminates the need to resolve the complicated angle 
dependencies within the integral of Eq. 6. Along with Eq. 7, Marshak's boundary 
conditions were used: 
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BG 

— ± 

Bz 


BC ? 
Br 

BG 

Br 

hG = 


0 ( Symmetry axis at r = r,) 

+ h G =4 n h l b ( wall at r = rj 

± 4 v;h/ b inlet Z=Z,(-) ex/fc Z=2^(+) 

Yi fc) 


( 8 ) 


For cases in which the boundary does not radiate back into the interior, the right- 
hand-side of the second and third lines in Eq. 8 become zero. Presently, the code 
sets the boundary emissivity € to 1 and the right-hand-side of Eq. 8 to 0 to 
represent complete absorption by the wall. However, provisions have been made 
to allow the boundaries to radiate back into the interior by setting a switch in the 
code, when e = 1, such that the RHS of Eq. 8 is evaluated. 

4.1.1 Radiation Coupling to FDNS 

The radiation source term was added to FDNS and calculated at every n 
global iterations, where n was chosen by computational experiments so as to 
reduce the computational expense of frequent coupling. Within each call to the 
radiation module, the current values of temperature, pressure and species number 
density obtained from the flow solution were used to generate new radiation 
properties (<7 ext , to, G) and then radiation intensity at each point. It is easier to 
solve for the radiation source term within physical space coordinates, after which it 
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can be re-cast in computational space coordinates quite simply. Also, radiation 
grid densities different from flowfield grid densities could be used by following this 
procedure. Following the procedure that transformed the energy equation from 
physical to computational space, and normalizing with reference values, the 
radiative source term is added to the right-hand side of the energy equation in 
FDNS as 


r- 1 


- v • Q = kT (4 n L - G) a A. 

Boltzmann Number * v b 


(9) 


where 


Boltzmann number = 


(pjjCpT^ 

Tt, 


° 'ref 

J is the transformation Jacobian, and o is the Stefan-Boltzmann constant. The 
Boltzmann number is the ratio of radiative to convective heating. No other terms 
need to be added to the Navier-Stokes equations since the following effects are 
negligible: radiation pressure, spectral doppler shift (due to material motion), and 
change in material thermal energy (due to material coupling). The source term is 
added to the interior grid points only; no special treatment is required for the 
heating rate boundary condition, as the gradients imposed at the wall points will 
adjust to the interior values. The subroutines that were developed for the ODA 
procedure are discussed in the Appendix. 


4.1.2 Gaseous Radiation Example 


Equations 7 to 9 were discretized and added to the FDNS code. A test case 
was run using a Rocketdyne motor which burned RP-1 and 0 2 . Once a converged 
solution was obtained without radiation, the radiation computations were initiated. 
No soot particles were considered in this case. The values of G were initialized 
with Planck functions at eleven wavelengths (H 2 0 and C0 2 ) at each flowfield grid 
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temperature. Eleven equations of the form of Eq. 7 were solved using an SOR 
method, with convergence occurring within 10 iterations. Once the value of G 
was computed at each interior point, the boundary values were determined 
according to Eq. 8 with boundary emissivity set to 1 and with no boundary re- 
radiation. The radiation computations were called every 50 iterations with a 0.1 
percent change in the source term from the previous call. 

Subsequently, the radiative source terms were evaluated and added to the 
energy equation. However, the solution did not show any appreciable difference in 
the temperature field due to the radiation source terms in either case. For the 
coupling radiation/convection nozzle flow solution, two limits were observed. In 
the high limit, the magnitude of the radiation source term was as large as ten times 
the non-radiative source term within the combustor, and in the low limit, as low as 
0.01 times the non-radiative source term (near the wall). In either situation, 
radiation still has a negligible effect on the flowfield. This may seem contradictory 
at first. However, this behavior follows from the way in which FDNS evaluates the 
source terms. The explicit or right-hand-side source terms include dissipation, 
conduction, the off-diagonal convection terms, wall function terms next to the 
wall, plus the radiation term. Diagonal convection terms are split off from the 
right-hand-side and moved to the left-hand-side to become implicit source terms. 
Since the implicit source terms in the energy equation dominate the explicit 
radiation term by a factor of 1e4 in the low limit, it is clear that the radiation 
source term is negligible compared to the other source terms. For cases where the 
implicit source term is zero, as in the high limit, the total explicit source term was 
negligible with or without radiation. 

The negligible radiation in the combustor is due to two effects: a constant 
temperature field and the large Boltzmann number of 7700. To test conditions for 
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which radiation could be important, a hypothetical test case was constructed: a 
linear variation in temperature was imposed across the motor, and the Boltzmann 
number was multiplied by a factor of 1e-4. Temperature contours of this start 
solution are presented in Fig. 5. After 2000 iterations of radiation coupling, the 
temperature contours take the form shown in Fig. 6. Comparison of radiation and 
non-radiation temperature contours indicates a reduction of temperature in the 
combustor and an increase of temperature in the nozzle with an average of five 
percent in each case, and that the source term was opposite in sign to that in the 
combustor due to the cold flow there. This test indicates that radiation can be 
expected to be important for cases where the temperature field differs from point 
to point, and where the convective heating rate does not completely dominate the 
radiation component. It can be concluded that radiation in the optically thick 
combustor region of the Rocketdyne test case was minimal at best. 

4.1.3 Particulate Radiation Example 

A coupled particulate radiation/CFD solution has not yet been attempted. All 
of the necessary Al 2 0 3 radiation subroutines have been added to FDNS, and 
computations can be made in the future to evaluate the effects of particles on 
flowfield properties. 

4.2 Improved Differential Approximation (IDA) 

The IDA was developed to extend the validity of the - approximation to 
also predict optically intermediate conditions. It involves a two-step process: the 
first is to obtain G from the P, - approximation, the second is to upgrade this 
solution by substitution into an additional partial differential equation (Ref. 2 and 
11). There is no integro-differential equation to solve; however, an evaluation of 
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several integrals is required. The following development will delineate the three 
steps required to obtain an IDA solution: i.) solution to ODA, ii.) solution to a 
surface integral equation for wall radiosities, and iii.) evaluation of surface 
integrals for points inside the medium. The source term to the ODA can be written 

S’(r.fi) = (1 -w) lb V) + t~ [G‘(r) + Pl q{r) ] (10) 

4ti 

Equation 10 was obtained from Eq. 4 after transforming the independent variables 
(r, z) to (y r , k z )- Also, the asterisk in Eq. 10 denotes the ODA or intermediate 
solution. 

A solution for \ A can be constructed from terms related to wall (w) and 
medium (m) contributions. 

I(r.&) ■ ur.fi) * /„(r,fi) ni) 

where 

T 

/ m (r,D) = f S*{r' + O'fi, (k' s H2) 

0 

and 

ntf) . MA e -'. (i3) 

n 

It is important to note that l w accounts for emission from the enclosure wall (that 
may have been attenuated by absorption and scattering in the medium and by 
reflections from the walls). I m is a function of the radiation released within the 
medium into a given direction by emission and scattering, and as such is related to 
the radiative source term (Ref. 2). 

The crux of the IDA is the way the source term is approximated. In optically 
thin regions, the scattering integral is negligible (the last two summations in Eq. 
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10) and attenuation follows a linear variation. These observations justify 
approximating the source term integral in Eq. 12 by a Taylor series. This results in 
the following representation for the medium: 

Ijr.ft) * S‘(r-o 0 a,fi) (i-e"'") ' 14) 

where t s refers to the optical distance between two points along a path length s (in 
this case from a point of interest in the medium to a point on the wall): 


x 


S 



( 15 ) 


The source term in Eq. 14 is evaluated at an adjusted point corresponding to an 
optical distance t 0 away from the point of interest, according to: 


T o ■ 1 - 




1-e 


( 16 ) 


This form for S’ is a consequence of using the assumed linearity of S’ in reverse, 
i.e., by linearly varying S" from the wall to the point under consideration. 


Radiosity from a particular wall point is composed of emission plus radiance 
on that point from other wall points, multiplied by the appropriate attenuation 
factor. Specifically, 

JJf) = € rt ! b jr) 


(i - £ ) / V. <rj • - * , um COS9COS9 » dA 


( 17 ) 


I A * nr ' W' w ' ' * * n o 

* nS 

The integral in Eq. 17 is evaluated as a summation over each wall segment j, 

where the view factors must be determined. A code or subroutine like the 

RAVFAC code (Ref. 12) may be used for this calculation. The attenuation factors 

involving e ™ are then evaluated as T tl between wall points i and j. Similarly, the 

source term takes the form S^ and corresponds to a value that lies between the 

wall points i and j according to Eq. 16. 
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The final step is to determine the improved value for the incident radiation G 
according to: 

G(r) = - / Wjrj e'- ♦ n ijrM a a < 18 > 

71 J ** 71 S 

The radiative source term is determined as in the ODA, using Eq. 9. Since this 
method is valid for all optical thicknesses, it can be used throughout the nozzle 
from combustor to exit, and in the plume. Most of the coding necessary to 
perform the IDA calculations has been made. A subroutine to implement the IDA 
has been written. This subroutine requires r s and view factors; methods to 
calculate these values will be developed in Phase II. A preliminary evaluation of r s 
and the required view factors was made by evaluating radiation in a cylindrical 
geometry configuration, as described in Section 4.3. 

4.3 Radiation Heat Transfer in an Infinite Cylinder 

Since view factors for an axisymmetric motor geometry were needed for 
implementing the IDA and since approximate methods of evaluating the optical 
depth between points in the motor, t s , are required, a cylindrical geometry was 
investigated. Rocket motor radiation was estimated with a diffusion approximation 
for a cylindrical geometry by Pearce (Ref. 13). Since LOS/view-factor analyses 
were to be made for a cylindrical geometry, the two-flux particle radiation model 
(Ref. 7) was used to provide an alternate solution to the problem Pearce posed. 
Such a solution involves determining the optical depth across the cylinder - data 
which can subsequently be used to judge how t s should be evaluated. The actual 
evaluation of r, will be a Phase II task. This approach also allows the development 
of a simple method for estimating rocket motor chamber wall radiation. 
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4.3.1 The LOS Method of Predicting Uncoupled Radiation 


An LOS method for predicting radiation heat transfer to a point on 
an infinite cylinder was developed by considering the cylinder to have constant 
temperature and constant composition. The total hemispherical radiation heat rate 
to a point on the edge of a cylinder can be found from the following equation: 


®rad 


fo ?fo* q < e ’*) d6d * 

/*/*cos0sIne dddfy 


(19) 


which can be reduced for the infinite cylinder to : 


rr^e.vdedt 

rad 


2n 


( 20 ) 


The term, q(0,0) , is the radiation heat flux to the point from a particular line 
of sight. It is evaluated in this line of sight code with a one-dimensional slab model 
using absorption and scattering OPTROCK data. Since the conditions inside the 
cylinder are considered isothermal and homogeneous, the radiation heat flux term 
for each different line of sight changes only with path length. Path lengths for a 
iine-of-sight inside an infinite cylinder are calculated as illustrated in Fig. 7, where 
Pi represents the path length for a individual line-of-sight. The heat rate terms for 
each different path length are numerically integrated according to Eq. 20 to obtain 
the total hemispherical radiation heat flux to the point. 


For Al 2 0 3 , particle radiation inputs into the program are the number and size 
of the particles, the particle density in the flow volume for each particle size, and 
temperature of each particle size. 
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4.3.2 Validation of the LOS Method 

To validate the program, a case was run to compare with Pearce (Ref. 13). 
The inputs which were used for this comparison are shown in Table 5. Particle 
sizes were chosen according to the particle size distribution used by Pearce. 

Particle densities were calculated by multiplying particle mass flow (calculated 
through the number density of each chosen particle size) by the density of the flow 
and dividing by the total mass flow. The total mass flow and density were 
calculated by running the CEC program for a rocket motor of similar propellant and 
chamber pressure. Table 6 lists the results of both methods for a cylinder with a 
radius of 3 inches. These results are in good agreement. Pearce also ran a series 
of cases with varying nozzle throat diameters assuming the same particle 
distributions. This data were compared to an empirical equation for determining 
emissivity in a solid rocket motor. The line-of-sight method, Pearce method and 
the empirical method results are all compared if Figure 8. The modest differences 
between the Pearce method and the line-of-sight method may be due to the 
particle size distribution differences, CEC rocket motor predictions differences, and 
different particle property data for the prediction of the absorption and scattering 
coefficients. 

Both the line-of-sight method and the Pearce method use the same 
assumptions of constant temperature and constant composition within the motor. 
Radiation predictions should be relatively accurate in the nozzle throat and 
combustion chamber where the assumed flow is accurate, as compared to the exit 
plane where temperature and density (composition) variations are more prominent. 

For supersonic portions of the nozzle, other methods which account for 
temperature and composition variations, as well as, geometrical considerations are 
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warranted. The LOS method is not restricted to constant flow properties or 
optically dense media. If a CFD grid and properties at the grid points are provided, 
an uncoupled LOS calculation through the temperature and composition gradients 
to the nozzle walls can be made. 


Table 5 

Input to Program for Comparison to Pearce [13] 

Particle No. 

Particle Size 
(Arm) 

Particle Mass/Flow Volume 
(slug/ft 3 ) x (10 4 ) 

Temperature 

(°R) 

1 

0.3 

2.8796 

5580. 

2 

0.5 

4.3929 

5580. 

3 

0.7 

4.7292 

5580. 

4 

1.0 

3.9305 

5580. 

5 

1.2 

3.0898 

5580. 

6 

1.5 

1.9758 

5580. 


Table 6 

Comparison of Predicted Radiation Heat Transfer @ R = 3 inches 

(BTU/ft 2 -s) 

Pearce [13] 

Line of Sight Program 

66.5 

70.7 
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5.0 Technical Feasibility 

This Phase I study established that a coupled radiation/convection CFD 
solution for rocket motor and plume analysis is feasible using variations of the 
differential approximation to describe the radiative transfer. The order and type of 
improvement required in the differential approximation depends on the specific 
application. Additional data storage is required to analyze radiation/convection 
coupled flowfields, but the additional iterations required to reach a solution with 
the modified FDNS code is modest. Further Phase II research is required to obtain 
the generality needed to analyze all types of propellant systems of interest and to 
provide a sufficient number of coupled solutions to determine when this 
methodology is preferable to conventional uncoupled methods. Specific 
conclusions drawn from the Phase I study are given by the following. 

The RTE was approximated by the P, - approximation method for application 
to optically thick media. A coupled radiation/convection solution was generated 
for an H 2 /0 2 rocket using the Pt - approximation to solve for the incident radiation 
at a point. This solution was driven by the extinction coefficients derived from a 
wide band model. Results indicate a negligible radiation effect in the combustor 
where the ODA is valid, for the one motor case considered, due to the very large 
convective effect indicated in the Boltzmann number. A hypothetical test case 
was run which indicated that the coupling calculation behaved properly. Particle 
radiation/flowfield coupling routines were completely coded, and additional runs 
will be made as part of a Phase II investigation. 

Extension of the ODA to the IDA was begun and is almost completed. The 
ODA solution is the first step in obtaining an IDA solution, which is valid 
throughout all optical thicknesses. In particular, the IDA will be used in the nozzle 
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region and in the combustor to compare the differences in the two approximations. 
Results are not anticipated to change in the optically thick combustor, but the 
nozzle section should show a different solution behavior. The IDA must be used 
for radiation/convection coupled plume analyses. 

The two-flux LOS method which was originally written for a one-dimensional 
plane-parallel slab geometry, was extended to cylindrical configurations. Optical 
thickness was computed for various angles across a single cylindrical slab. Good 
comparisons for radiation heating predictions were obtained between LOS and a 
diffusion approximation for the nozzle throat region. For optically thin regions, 
more slabs are required but the LOS method is still applicable. Furthermore, the 
LOS method can be used to advance and evaluate the development of the IDA 
method. 

Finally, the radiation updates made to the FDNS code to provide the 
radiation/convection coupled solution are presented as modules on a floppy disk 
delivered to the technical monitor of this work. A description of the radiation 
module updates to FDNS is given in the appendix. 
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Appendix 

A brief description of the radiation updates added to FDNS is given in this 
section. This code will be provided on a floppy disk supplied to the technical 
monitor of this work along with the final report. 

Five of the original FDNS subroutines were upgraded to include radiation call 
statements and include files. Their names were appended with _RAD to signify 
these changes. 


DATAIO RAD: 


include fdns20 
include dataio3.inc 
include dataio4.inc 


FDNS RAD: 


include fdns20 
include fdns3.inc 
include fdns4.inc 
call radO 


PSOURC RAD: 


include fdns20 

eliminates radiation emission term of original 
module since its effect is added through the 
radiation source term in radinO. 


SOURCE RAD: 


include fdns20 
include source. inc 


I NIT RAD: 


include fdns20 

include init6. inc 

reads in radiation switch IRAD 


New Radiation Modules 


Fortran Files 

astar dimensionless band absorption for H 2 0 and C0 2 
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Fortran Fi 

ecoefg 

ecoefp 

extinc 

htrad 

linint 

logint 

pathl 

plankf 

ptlrad 

radO 

radbc 

radinO 

radpl 

refrac 

lininr 

refrin 

sigacl 

sigav 


(Cont.) 

extinction coefficient for H 2 0 and C0 2 

extinction coefficient for Al 2 0 3 

driver for extinction coefficients 

radial and axial heating rate at a point (Btu/ft 2 -s) 

double linear interpolation of particle size and temperature for the 
absorption coefficient 

double log (base 10) interpolation of particle size and temperature for 
the scattering coefficient 

path length at each grid location 

Planck function evaluation 

transforms FDNS particle flowfield results to particle temperature and 
number density 

driver for radiation module 

radiation boundary conditions 

driver for P, - approximation; computes radiation source term 
SOR routine for incident radiation G 

reads in particle refractive index (used in evaluation of Planck 
function) from KOCH2.PRN data 

linear interpolation of refractive index as function of temperature 
driver for refractive index interpolation 

particle absorption coefficients from OPTROCK data located in 
SIGACL.DAT 

particle absorption and scattering coefficients 
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Fortran Files (Corn.) 

sigscl particle scattering coefficients from OPTROCK data located in 

SIGSCL.DAT 


summal and 

summa2 wide band model summation - function of vibrational quantum number 


trdiag tridiagonal solver for P r approximation 

widebm wide band model for H 2 0 and C0 2 
Data Files 


BLKDAT.F 

KOCH2.PRN 

SIGACL.DAT 

SIGSCL.DAT 


radiation data for H 2 0 and C0 2 

optical property data from Grumman OPTROCK test for Al 2 0 3 

absorption coefficients for Al 2 0 3 obtained from transforming 
KOCH2.PRN by Mie code 

scattering coefficients for Al 2 0 3 obtained from transforming 
K0CH2.PRN by Mie code 


Include Files 

fdns20 

fdns21 

fdns22 


alumox.inc 

h2o.inc 

co2.inc 

source. inc 

dataio3.inc 

dataio4.inc 

init6.inc 

fdns3.inc 

fdns4.inc 
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Initialization Modules 
Fortran File 

grid initializes grid and flowfield for nozzle flow without radiation 

Include Files 
grid 1 .inc 

grid2.inc initializes normalized incident radiation per wavelength and the 
refractive index; writes to units 15 and 16 respectively 

grid3.inc refractive index read and interpolation modules 

no provision has been made to read in separate flowfield and radiation 
grid files 

Figure 9 presents a flow chart of the radiation driver, subroutine radO.f. 
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